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TESTS  OF  PARAMETERIZED  LANGMUIR-CIRCULATION  MIXING  IN 
THE  OCEAN’S  SURFACE  MIXED  LAYER  II 


1.  INTRODUCTION 

It  has  long  been  assumed  that  Langmuir  circulation  (LC)  affects  mixing  in  the  ocean’s  sur¬ 
face  mixed  layer  (SML),  but  until  computing  power  had  advanced  sufficiently  to  allow  large-eddy 
simulations  (LES)  of  LC,  there  was  not  much  quantitative  information  upon  which  to  base  a  pa¬ 
rameterization  of  the  effects  of  LC  mixing  in  ocean  models. 

However,  beginning  in  the  1990’s,  LES  began  being  used  to  investigate  LC  and  its  effects 
on  upper-ocean  mixing.  Some  papers  on  this  research  include  Skyllingstad  and  Denbo  (1995), 
McWilliams  et  al.  (1997,  hereafter  referred  to  as  MW97),  Skyllingstad  et  al.  (1999),  Skyllingstad 
(2000),  McWilliams  and  Sullivan  (2000),  Min  and  Noh  (2004),  Noh  et  al.  (2004),  Sullivan  et  al. 
(2004),  Harcourt  and  D’Asaro  (2007),  Tejada-Martinez  and  Grosch  (2007),  Polton  and  Belcher 
(2007),  Sullivan  et  al.  (2007),  Grant  and  Belcher  (2009),  Kukulka  et  al.  (2009),  Alan  et  al.  (2009), 
Li  et  al.  (2009),  Tejada-Martinez  et  al.  (2009),  Martinat  et  al.  (2011),  Sullivan  et  al.  (2012),  Van 
Roekel  et  al.  (2012),  and  Hamlington  et  al.  (2014). 

MW97  conducted  LES  simulations  of  LCs  for  a  simple,  wind- mixing  case  with  an  initial  mixed- 
layer  depth  (MLD)  of  about  33  m  and  a  moderate  wind  of  about  5  m/s.  They  found  that  the 
inclusion  of  LCs  in  their  LES  increased  the  rate  of  mixing  within  the  SML  by  about  a  factor  of  3 
and  slightly  increased  the  depth  of  mixing. 

Based  on  these  results  from  MW97,  Kantha  and  Clayson  (2004,  hereafter  referred  to  as  KC04) 
investigated  how  to  implement  the  increased  rate  of  mixing  by  the  LC  in  the  widely-used  Mellor- 
Yarnada  Level  2.5  (MYL2.5)  turbulence  model  (Mellor  1974,  Mellor  and  Yarnada  1982).  They 
added  additional  shear-production  terms  to  the  MYL2.5  model,  consisting  of  the  product  of  the 
vertical  Reynolds  stress  multiplied  by  the  vertical  shear  of  the  Stokes  drift  current  (SDC)  from  the 
waves.  The  additional  shear-production  terms  were  added  to  both  the  turbulent  kinetic  energy 
(TKE)  equation  and  the  turbulence  length-scale  (TLS)  equation  of  the  MYL2.5  model.  For  the 
conditions  of  the  simple,  wind-mixing  simulation  conducted  by  MW97,  KC04’s  modifications  sig¬ 
nificantly  increased  the  maximum  rates  of  mixing  in  the  SML,  and  the  higher  mixing  rates  were 
fairly  consistent  with  the  LES  results  of  MW97. 

Martin  et  al.  (2013,  hereafter  referred  to  as  M13)  described  the  implemetation  of  KC04’s 
parameterization  of  LC  mixing  in  the  version  of  the  MYL2.5  turbulence  model  that  is  used  in  the 
Navy  Coastal  Ocean  Model  (NCOM)  and  the  results  of  a  number  of  tests  of  the  implementation  of 
the  KC04  parameterization  in  NCOM.  The  tests  included  the  simple,  wind-mixing  test  case  used 
by  MW97  and  KC04,  a  year-long  simulation  of  the  SML  at  Ocean  Weathership  Station  (OWS) 
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Papa  in  the  NE  Pacific,  and  a  simulation  of  the  effect  of  Hurricane  Ivan  on  the  upper  ocean  in 
the  Gulf  of  Mexico  in  2004.  The  results  were  found  to  be  consistent  with  those  reported  by  KC04 
in  that  their  parameterization  of  LC  mixing  in  the  MYL2.5  turbulence  model  increased  the  rate 
of  mixing  within  the  SML  by  a  factor  of  2-3  and  slightly  increased  the  depth  of  mixing,  and  was 
otherwise  well  behaved  and  numerically  robust. 

Since  the  KC04  paper,  there  have  been  additional  efforts  to  parameterize  the  effects  of  LC  on 
upper-ocean  mixing  in  turbulence  models.  Two  notable  efforts  are  those  by  Harcourt  (2013  and 
2015).  Harcourt  noted  that  the  stability  functions  used  in  KC04  were  derived  from  the  algebraic 
Reynolds  stress  model  (ARSM)  used  in  Kantha  and  Clayson  (1994,  hereafter  referred  to  as  KC94) 
that  included  the  effects  of  local  stratification  and  shear,  but  not  the  Craik-Leibovitch  (CL)  vortex 
force  (Craik  and  Leibovich  1976),  and  that  this  was  inconsistent  with  KC04’s  use  of  the  CL  vortex 
force  to  derive  modifications  to  the  prognostic  equations  for  predicting  the  TKE  and  the  TLS. 
Hence,  Harcourt  included  the  CL  vortex  force  in  all  the  ARSM  equations  used  to  determine  the 
effect  of  LC  on  both  the  stability  functions  and  the  TKE  and  TLS  equations. 

In  this  report,  the  implementations  in  NCOM’s  MYL2.5  turbulence  model  of  the  parameteri¬ 
zation  of  LC  mixing  by  Kantha  and  Clayson  (2004)  and  by  Harcourt  (2015,  hereafter  referred  to  as 
H15)  are  described,  and  the  results  of  some  tests  are  reported,  including  comparisons  with  NCOM’s 
original  MYL2.5  turbulence  model  without  a  parameterization  of  LC  mixing.  Hence,  this  report  is 
a  follow  on  to  the  discussion  of  the  KC04  parameterization  of  LC  mixing  in  M13. 

The  following  sections  contain  a  description  of  the  implementation  of  both  the  KC04  and  the 
H15  Langmuir  Circulation  mixing  parameterizations  (LCMPs)  in  the  MYL2.5  turbulence  model 
used  in  NCOM  (Section  2),  a  test  of  the  KC04  and  H15  LCMPs  for  the  simple,  wind-mixing,  test 
case  used  by  MW97  and  KC04  (Section  3),  a  test  of  the  KC04  and  H15  LCMPs  at  OWS  Papa 
(Section  4),  a  test  of  the  KC04  and  H15  LCMPs  for  a  simulation  of  Hurricane  Ivan  in  the  Gulf  of 
Mexico  (Section  5),  and  a  summary  (Section  6). 

2.  IMPLEMENTATION  OF  LC  MIXING  IN  THE  MYL2.5  MLM 


2.1  Original  MYL2.5  MLM  in  NCOM 


The  original  TKE  and  TLS  equations  as  implemented  in  the  MYL2.5  MLM  used  in  NCOM 
are 


dq 2 
~dt 


-V  •  (yq2)  +  Qq 2  +  Vh(AHVhq2)  +  J- 


+  2  Km 


+  2  Kh 


g_dp_ 
Po  dz 


dq2£ 

1 W 


-V  •  ( vq2£ )  +  Qq2£  +  Vh(AHVh(q2£))  +  |- 


+  E1£Km 


+  E3lKH<Lf 

Po  OZ 


Q 

E2j-W. , 
b  i 


(1) 


(2) 


where  q  is  the  square  root  of  twice  the  TKE,  £  is  the  turbulent  length  scale,  t  is  the  time,  v  is 
the  vector  velocity,  Q  is  a  volume  flux  source  term,  Ah  is  the  horizontal  mixing  coefficient,  z  is 
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the  vertical  coordinate,  Kq  is  the  vertical  diffusion  coefficient  for  q 2  and  q2£,  which  is  taken  to  be 
proportional  to  Km  (Kq  =  OAIKm),  Km  is  the  vertical  mixing  coefficient  for  momentum,  u  and 
v  are  the  horizontal  components  of  the  velocity,  Ku  is  the  vertical  diffusion  coefficient  for  scalar 
fields,  g  is  the  acceleration  of  gravity,  pQ  is  a  reference  density  for  the  water,  and  dp/dz  is  the 
vertical  gradient  of  the  water  density,  with  the  contribution  of  the  pressure  to  the  vertical  density 
gradient  removed. 


W  is  referred  to  by  Mellor  and  Yarnada  (1982)  as  a  “wall  proximity”  function,  which  is  used 
to  scale  £  near  the  surface  and  bottom.  This  function  is  defined  by 


W  =  1  +  E4 


(3) 


and  in  NCOM,  L  is  defined  by 


L  1  —  (C  —  z  +  zs)  1  +  (z  —  H  +  Zb)  1, 


(4) 


where  k  =  0.4  is  Von  Karman’s  constant,  £  is  the  elevation  of  the  free  surface,  zs  is  the  surface 
roughness,  Zb  is  the  bottom  roughness,  and  H  is  the  static  bottom  depth.  The  symbols  bi,  E\ ,  E2 , 
E3,  and  E4  are  constants,  with  values  specified  in  Table  1. 

Table  1  —  Constants  for  MYL2.5  Turbulence  Equations 
used  in  NCOM 


parameter 

value 

ai 

0.92 

bi 

16.6 

CL2 

0.74 

b2 

10.1 

Cl 

0.08 

Ei 

1.8 

e2 

1.0 

E3 

1.8 

Ea 

1.33 

The  vertical  mixing  coefficients  Km  and  Ku  are  computed  as 

Km  =  £qSM, 

(5) 

Kh  =  £qSn, 

(6) 

where  Sm  and  Sh  are  “stability  functions”  computed  as 

Ci 

H  1  -C2Gh' 

(7) 

0  C3  +  CaGhSh 

Sm~  1  -  C5Gh  ’ 

(8) 
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where 


Gh  =  min 


0.028, 


£2g  dp 


(9) 


q2p0  Sz 

The  constants  C\  -  C$  are  calculated  from  the  basic  turbulence  constants  (ai,  02,  61,  b2,  and  c\)  as 

Ci  =  a2(bi  -  6ai)/bi,  (10) 


C*2  —  a2(18ai  +  362)1 


(11) 


C3  =  ai(bi(l  -  3ci)  -  6ai)/6i, 


(12) 


C4  =  ai(18oi  +  902), 


(13) 


C5  —  90102. 

Table  1  lists  the  values  of  the  constants  used  in  the  MYL2.5  turbulence  model  in  NCOM. 


(14) 


2.2  Implementation  of  KC04  parameterization  of  LC  mixing  in  MYL2.5  MLM 

KC04  parameterized  the  effects  of  the  LC  mixing  from  the  LES  conducted  by  MW97  by  adding 
additional  shear  production  terms  to  both  the  TKE  and  TLS  equations  in  the  MYL2.5  turbulence 
model.  The  TKE  and  TLS  equations  in  NCOM  with  the  additional  shear  production  terms  used 
by  KC04  are 


dq2 

~dt 


=  -v  •  (W)  +  Qq2  +  vh(AHvhq2)  + 


d_ 

dz 


K, 


dq 2 
dz 


+  2  Km 


du\  2  ( dv\2  du  dus  dv  dvs  \  „  q  dp  „  g3 

at)  +(&)  +Tz ~at  +  lTz~at  )+2A"y&“2w’ 


(15) 


dq2i 

~df 


d 


dz  \  “q  dz 


-V  •  (vq2£)  +  Qq2£  +  Vh{AHVh(q2l))  +  £  (  K, 
+  EiIKm 


dq2£ 


du\2  { dv\2\  „  fdudus  dv  dvs 

at)  +  (at)  )  +  E“eKM \at~Bz  +  Tz~at 

9  dp  „  q3 


+  E3£Kh^^-E2^-W, 

Po  dz  b ! 


(16) 


where  us  and  vs  are  the  horizontal  components  of  the  SDC  from  the  waves.  An  additional  constant 
Eq  is  used  to  scale  the  Stokes  shear  production  term  in  the  TLS  equation  and  has  a  value  of  7.2. 
Note  that  the  value  of  7.2  for  Eq  is  different  from  the  value  of  4.0  reported  in  KC04,  which  was 
incorrect  (Dr.  Lakshmi  Kantha,  personal  communication).  We  found  in  early  testing  that  the  KC04 
LC  mixing  parameterization  does  not  work  correctly  (i.e.,  as  reported  by  KC04)  with  E$  =  4.0. 

The  original  shear-production  terms  are  proportional  to  the  square  of  the  vertical  shear  of 
the  Eulerian  velocity  u  and  v,  hence,  the  original  shear-production  terms  are  always  positive  and 
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always  act  to  increase  the  TKE.  The  shear-production  terms  added  by  KC04  to  parameterize  the 
LC  mixing  include  the  vertical  shear  of  the  model  Eulerian  velocity  times  the  vertical  shear  of  the 
SDC;  hence,  it  is  possible  for  this  term  to  be  negative  if  the  two  velocity  shears  are  of  different 
sign,  which  would  result  in  a  loss  of  TKE  from  this  term,  rather  than  a  gain.  In  tests  conducted  to 
date,  this  term  is  usually  positive  when  the  winds  and  wave  directions  are  approximately  alligned. 
When  the  winds  and  waves  are  in  significantly  different  directions,  this  term  can  be  negative  and, 
hence,  act  to  reduce  the  amount  of  mixing.  A  reduction  in  the  vertical  mixing  with  increasing 
misallignment  of  the  wind  and  wave  directions  has  been  shown  to  occur  in  LES  (Van  Roekel  et 
al.  2012).  Preliminary  tests  with  the  KC04  (e.g.,  Carniel  et  al.  2005)  and  H15  LCMPs  have  also 
shown  a  reduction  in  vertical  mixing  with  increasing  misallignment  of  the  wind  and  waves,  though 
determination  of  the  extent  of  the  agreement  of  the  LCMPs  with  the  LES  in  this  regard  needs 
further  investigation. 

The  vertical  mixing  coefficients  and  stability  functions  are  computed  using  the  equations  de¬ 
scribed  in  the  previous  section.  The  constants  used  by  KC04  are  the  same  as  those  listed  in  Table  1 
in  the  previous  section,  except  that  E4  =  4.87. 

2.3  Implementation  of  H15  parameterization  of  LC  mixing  in  MYL2.5  MLM 


The  derivation  of  H15’s  parameterization  of  the  effects  of  LC  on  upper-ocean  mixing  can 
be  found  in  Harcourt  (2015).  What  will  be  described  here  is  the  implementation  of  the  H15 
parameterization  of  LC  mixing  in  NCOM. 


H15’s  inclusion  of  the  CL  vortex  force  in  the  ARSM  equations  results  in  a  modified  form  of  the 
vertical  turbulent  momentum  flux,  i.e., 


K  9U  _1_  K*  9Us 

AM^  +  AM  — , 

(17) 

K  dv  4.  T<S  dVs 

KmTz  +  Km  aP 

(18) 

where  the  second  term  on  the  right-hand  side  (RHS)  of  these  two  equations  represents  a  vertical 
turbulent  momentum  flux  down  the  gradient  of  the  Stoke’s  current  proportional  to  the  vertical 
mixing  coefficient  Kfr.  This  additional  term,  which  does  not  appear  in  the  KC04  parameterization, 
affects  the  momentum  equations  for  u  and  v,  as  well  as  the  shear  production  terms  in  the  TKE  and 
TLS  equations.  When  the  momentum  equations  are  solved,  this  additional  term  acts  as  an  explicit 
forcing  term  on  the  RHS  of  the  momentum  equations. 


The  TKE  and  TLS  equations  with  the  additional  shear  production  terms  derived  by  H15, 
including  the  modified  vertical  momentum  flux,  are 


dq 2 
dt 


d 


dq 2 


-V  •  (V)  +  Qq2  +  Vh(AHVhq2)  +  —  (  Kq-A- 


+  2 


du 

dz 


du 


+  Hp0dz  2m 


M  dz 
3 


du 

dz 


km  —  +  ~  +  -tt"  )  +  (  km—  +  ^m~cT~ 


dz 

dus 

dz 


dv 

dz 


dv. 

dz 


/<9u  dvs\\ 
\dz  +  dz  JJ 


(19) 
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dq2 


dt 


=  -V  •  (vq2£)  +  Qq2i  +  Vh(AHVh{q2£ ))  + 


d_ 

dz 


K 


dq 2 


8z 


<9it  chts\  <9«  /  .  <9u  ch/O  <9f 

+  I  (  KmTz  +  ^  +  ks  +  Am^J  ^ 


<9tt 


9u,\  <9« 


cO 


ch/,\  cO 


+  Ee£  (  (  Km t—  +  Em— A-  )  — -  +  (  KM~z~  +  K%—^~  )  -— 4 


9^  '  dz  7  32  l  M  dz  J  dz 

g  dp  ^  r 


+  E3IKh-^  -E2j-W. 

Po  dz  b  1 


(20) 


The  vertical  mixing  coefficients  are  computed  as 

Km  =  £qSm, 

Kh  =  £qSn, 

Km  =  ^qSm/z 


(21) 

(22) 

(23) 


Kq  =  0.41/0/, 


(24) 


where  fz  is  a  “surface  proximity  function”,  which  is  an  empirical  function  chosen  by  H15  to  reduce 
values  near  the  surface  in  order  to  improve  agreement  with  observations. 


The  stability  functions  are  computed  as 


ci>  _ 
z>m  ~ 


C  i 


1  -  C2Gh  -  C:iGvfz 


(25) 


Sh  = 


CO  -  C12Gh  +  C13Gsf 2  -  CuGvfz 


(1  -  C15Gh  -  C16Gvfz)(  1  -  C17Gh)  -  (Ci8  +  C19GHfz  -  C29Gvfz)Gvfz 


(26) 


where 


SM  = 


C31  +  C32GhSh  +  C33Gsf2SsM 


1  -  CUGH  -  C35Gvfz 


KqJ  p0dz 


(27) 


(28) 


Gv  = 


£  \  ^  du  dus  +  dv  dv< 


qj  \dz  dz  dz  dz  J  ' 


(29) 


2  /  '  duR  \  2  .  fdvs  V 


dz 


+ 


dz  J 


(30) 
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The  surface  proximity  function  fz  is  computed  as 

fz  =  tanh  {Cfz/is), 


where 


,S  fJPldz 

L  H**  ’ 


where  Cf  =  0.25,  Pf  =  max(0,P),  and  P  is  the  Stokes  shear  production  term 


-dUe 


-dv« 


P  =  — u'w - v'w'  -  „  , 

dz  dz 


which  is  computed  as 


du  „s  dus\  dus 


dz 


P  =  ^  SM-  +  S^  SM-  +  S^  ^  . 


dz  )  dz 


dv  „s  dvs\  dvs 


dz 


dz  )  dz 


(31) 


(32) 


(33) 


(34) 


The  surface  proximity  function  is  small  (but  positive)  near  the  surface  and  increases  to  a  value  of 
one  within  the  SML. 


Because  of  the  dependence  of  the  surface  proximity  function  fz  on  the  stability  functions,  the 
calculation  of  fz  and  the  stability  functions  is  iterated  a  few  times  (5  times  in  the  simulations 
conducted  for  this  report)  to  obtain  approximate  convergence.  The  order  of  these  calculations  in 
NCOM  is  roughly  the  inverse  of  the  order  of  the  equations  above,  i.e.,  (34),  (32),  (31),  (28-30),  and 
(25-27).  After  the  stability  functions  are  computed,  new  values  of  the  vertical  mixing  coefficients 
(21-23)  can  be  calculated. 

The  form  of  the  stability  functions  presented  in  H15  involves  more  constants  than  those  pre¬ 
sented  in  (25-28).  Equations  (25-28)  are  the  calculations  used  in  NCOM  and  represent  the  maxi¬ 
mum  reduction  of  H15’s  stability  functions  that  can  be  obtained  by  combining  constants.  Relations 
between  the  constants  used  in  H15  and  those  in  (25-28)  are  presented  in  Table  3. 

The  complexity  of  H15’s  stability  functions  allows  for  the  possibility  of  unreasonable  or  unstable 
behavior,  e.g.,  the  denominators  of  (25-28)  could  possibly  go  to  zero.  Hence,  as  noted  in  H15,  some 
“realizability  constraints”  are  needed  to  maintain  reasonable  values  and  avoid  numerical  instability. 
Since  no  guidance  regarding  realizability  constraints  were  included  in  H15,  the  following  constraints 
were  implemented  in  NCOM  based  on  some  trial  calculations:  the  maximum  value  of  Gh  was 
limited  to  0.032,  the  denominators  of  (25-28)  were  restricted  to  minimum  values  of  0.01,  and  the 
stability  functions  themselves  were  limited  to  values  in  the  range  of  0  to  5.  Also,  the  vertical  mixing 
coefficients  (25-27)  were  limited  to  the  range  of  0-10  m2/s.  Note  that  not  a  lot  of  time  was  spent 
investigating  this;  hence,  these  constraints  could  no  doubt  be  improved. 

The  values  of  the  constants  used  by  H15  are  listed  in  Tables  2  and  3.  Note  that  E%  =  5.0 
instead  of  1.8  as  used  in  KC04  and  Eq  =  6.0  instead  of  7.2. 

The  behavior  of  the  H15  turbulence  parameterization  when  the  SDC  is  zero  is  of  interest,  since 
the  SDC  in  the  SML  will  be  small  when  the  waves  are  small,  and  will  be  zero  in  deep  mixed  layers 
below  the  influence  of  the  surface  waves.  When  the  SDC  is  zero,  the  Gy  and  Gs  functions  ((29) 
and  (30))  will  be  zero,  and  the  stability  functions  Sh  and  Sm  become 
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Table  2  —  Constants  for  Turbulence  Equations  for  H15 


parameter 

value 

a\ 

0.92 

h 

16.6 

C12 

0.74 

b2 

10.1 

Cl 

0.08 

C2 

0.7 

C3 

0.2 

Ei 

1.8 

E2 

1.0 

e3 

5.0 

Eii 

1.33 

e6 

6.0 

Table  3  —  Constants  for  Stability  Functions  for  H15 


parameter 

value 

equivalent  base  turbulence  constants 

7i 

0.667470 

1  —  607/61 

Ci 

0.614072 

ai7i 

c2 

6.127200 

9aia2 

c3 

7.617600 

9a  j 

Cn 

0.493928 

0271 

Ci2 

3.026394 

9aia27i 

Cl3 

10.551480 

9a1a271(2a1  +  a2) 

Cm 

8.292656 

9aia2(2ai(7i  +  3ci)  -  a2(7i  -  3ci)) 

Cm 

6.127200 

9aia2 

Cl  6 

30.470400 

36  aj 

C17 

30.192001 

3a2(6ai  +  62(1  -  c3)) 

Ci8 

1.478520 

9d2(l  -  c2) 

Cig 

284.308960 

162afa2(2ai  +  a2(2  -  c2)) 

C20 

45.051102 

324afa2(l  —  c2) 

C31 

0.393272 

u 

CO 

1 

e 

C32 

17.073360 

9ai(2ai  +  a2(l  —  c2)) 

C33 

22.852800 

27  af 

C34 

6.127200 

9aia2 

C35 

30.470400 

36  aj 

Val  ue 
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Stability  Functions  SH  Stability  Functions  SM 


Fig.  1  Plots  of  stability  functions  for  MYL2.5  and  H15  with  SDC  set  to  zero. 


Sh 


C 11  -  C12GH 
(l-C15GH)(l-C17GHy 


(35) 


Sm  = 


C31  +  C32GhSh 
1  ~  CuGh  ' 


These  can  be  compared  with  the  stability  functions  for  the  MYL2.5  MLM  in  (7)  and  (8).  The 
equations  for  Sm,  (8)  and  (36),  have  the  same  form  and  the  constants  are  the  same  except  for 
C4  for  MYL2.5  and  C32  for  H15,  which  though  different  are  somewhat  similar  in  value,  i.e. ,  21.36 
and  17.07,  respectively.  The  equations  for  Sh,  (7)  and  (35)  have  slightly  different  forms,  but  the 
constants  are  such  that  the  function  values  are  similar.  This  is  illustrated  in  Fig.  1,  which  shows 
plots  of  Sh  and  Sm  for  MYL2.5  and  for  H15  with  the  SDC  set  equal  to  zero.  Hence,  it  might  be 
expected  that  when  the  SDC  is  zero,  H15  will  predict  vertical  mixing  similar  to  MYL2.5;  this  is 
demonstrated  in  Section  4.4  in  simulations  at  OWS  Papa. 
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2.4  Additional  changes  to  NCOM  to  accommodate  wave  forcing 

The  coupling  of  NCOM  with  the  output  from  a  wave  model  requires  some  additional  changes 
besides  enhancing  the  vertical  mixing  within  the  MYL2.5  turbulence  model.  The  SDC  is  added  to 
NCOM’s  Coriolis  term  as  described  by  MW97  and  KC04  and  many  others.  This  term  is  referred  to 
as  the  Stokes-Coriolis  term  and  it  affects  the  downwind  transport  in  the  surface  wind-driven  layer 
as  discussed  in  MW97  and  KC04  and  in  this  report  in  the  next  section.  The  SDC  is  also  used  to 
advect  the  model  fields  in  NCOM  and,  to  be  consistent  with  this  advection,  is  added  to  NCOM’s 
continuity  equation.  Note  that  the  addition  of  the  SDC  to  NCOM’s  Coriolis  terms  affects  local, 
one-dimensional  simulations,  such  as  those  presented  in  this  report,  but  the  advection  by  the  SDC 
does  not. 

Additional  wave  effects  include  the  surface  wave-radiation  stress,  which  is  implemented  in 
NCOM  similar  to  the  surface  wind  stress;  but  since  it  is  generally  much  smaller  than  the  wind  stress, 
its  effect  is  usually  small,  though  not  negligible.  Another  effect  of  the  waves  is  the  enhancement  of 
bottom  drag  in  shallow  water  due  to  the  orbital  motions  of  the  surface- waves  near  the  bottom  (Grant 
and  Madsen  1979;  Soulsby  1995).  Both  of  these  additional  wave- forcing  effects  are  implemented  in 
NCOM,  but  these  effects  are  not  the  focus  of  this  study  and  will  not  be  discussed  further  in  this 
report. 

3.  TEST  OF  LC  MIXING  IN  MYL2.5  MLM  FOR  SIMPLE  WIND-MIXING  CASE 

MW97  used  a  simple,  wind-mixing  test  case  to  look  at  the  effects  of  LC  on  vertical  mixing,  and 
KC04  ran  the  same  case  to  test  their  parameterization  of  LC  mixing  in  the  MYL2.5  MLM  and  to 
compare  their  results  with  those  of  MW97.  Hence,  we  ran  the  same  case  to  test  our  implementations 
of  KC04’s  and  H15’s  parameterizations  of  LC  mixing  in  the  version  of  the  MYL2.5  MLM  used  in 
NCOM  and  to  compare  with  the  results  of  MW97  and  KC04. 

The  initial  condition  for  this  test  case  consists  of  a  SML  with  a  depth  of  33  m  and  a  temperature 
of  13.5°C  and  a  thermal  stratification  below  the  SML  of  0.01°C/m.  The  salinity  is  set  to  a  uniform 
35  psu  both  within  and  below  the  SML.  The  initial  temperature  was  chosen  to  give  a  thermal 
expansion  coefficient  of  -0.0002  1/°C.  The  latitude  is  43.4°N  to  give  a  value  of  the  Coriolis  parameter 
of  0.0001  1/s,  which  results  in  an  inertial  period  of  17.45  h. 

The  forcing  for  the  test  case  consists  of  a  surface  wind  stress  of  0.037  Pa,  which  corresponds  to 
a  wind  speed  of  about  5  m/s.  MW97  used  a  small  cooling  surface  heat  flux  of  5  W/m2  to  speed  up 
the  convergence  of  their  LES  numerical  solution  and  reduce  its  run  time.  We  used  a  linear  ramp 
the  length  of  the  inertial  period  (17.45  h)  for  the  surface  wind  stress  to  reduce  inertial  oscillations 
and  speed  the  convergence  to  an  approximately  steady  solution. 

The  surface  wave  field  for  this  test  case  was  specified  by  MW97  as  a  monochromatic  wave  of 
amplitude  0.8  m  and  wavelength  60  m  propagating  in  the  downwind  direction.  This  gives  a  SDC 
with  a  surface  value  of  0.0679  m/s  and  an  e-folding  depth  (the  inverse  of  twice  the  wavenumber, 
see  Appendix  A)  of  4.78  m,  which  corresponds  to  a  light  swell. 

NCOM  was  run  for  this  case  in  what  we  refer  to  as  pseudo  one-dimensional  mode  by  running 
the  model  on  a  horizontal  domain  of  2  by  2  grid  points  (the  minimum  horizontal  grid  that  NCOM 
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TKE 


KM 


U 


TKE  (cm2/s2) 


KM  (cm2/s)  U  (cm/s) 


V 


V  (cm/s) 


Fig.  2  —  Profiles  of  TKE,  vertical  mixing  rate  for  momentum,  and  downwind  and  crosswind  velocity  for  simulations 
without  a  LCMP  (solid  line)  and  with  the  KC04  LCMP  (dotted  line)  and  H15  LCMP  (dashed  line)  for  the  MW97, 
wind-mixing,  test  case. 


allows)  with  doubly-periodic  lateral  boundary  conditions.  The  vertical  grid  consisted  of  40  layers 
with  an  upper-layer  thickness  of  1  m  and  a  uniform  stretching  to  a  maximum  depth  of  200  m; 
hence,  each  layer  is  thicker  than  the  layer  above  it  by  about  7%. 

With  the  use  of  a  linear  ramp  for  the  surface  wind  stress,  the  numerical  solution  is  fairly  steady 
once  the  ramp  ends  after  17.45  h.  Figure  2  shows  profiles  of  TKE,  the  vertical  mixing  coefficient 
for  momentum,  and  the  downwind  and  crosswind  velocity  for  simulations  without  a  LCMP  and 
with  the  LCMPs  of  KC04  and  H15.  This  figure  can  be  compared  with  Fig.  2,  3,  and  4  in  MW97 
and  with  Fig.  1  in  KC04. 

The  TKE  in  Fig.  2  is  roughly  doubled  by  the  use  of  the  KC04  LCMP,  which  is  qualitatively 
similar  to  the  result  in  KC04.  For  the  H15  LCMP,  the  TKE  is  roughly  tripled,  except  for  the 
reduction  of  the  TKE  for  H15  in  the  upper  3  m.  Note  that  the  increase  in  TKE  in  the  upper  meter 
in  Fig.  2  for  the  original  MYL2.5  MLM  and  for  the  KC04  LCMP  is  due  to  the  parameterization 
of  the  surface  flux  of  TKE  from  the  surface  waves.  The  reduction  of  the  TKE  near  the  surface  for 
the  H15  LCMP  is  due  to  the  effect  of  H15’s  surface  proximity  function. 

For  the  KC04  LCMP,  the  maximum  value  of  the  vertical  mixing  coefficient  for  momentum  Km 
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in  Fig.  2  is  roughly  doubled  relative  to  the  original  MYL2.5  MLM,  which  is  qualitatively  similar 
to  the  result  reported  in  KC04,  i.e. ,  the  maximum  value  of  Km  for  the  KC04  LCMP  in  Fig.  2 
increases  from  about  200  to  400  cm2/s,  compared  with  an  increase  from  about  215  to  430  cm2/s 
in  KC04  (their  Fig.  1).  For  MW97  (their  Fig.  3b),  the  maximum  value  of  Km  increases  by  about 
a  factor  of  3  from  about  130  to  380  cm2/s.  However,  MW97’s  maximum  value  with  the  LC  mixing 
(380  cm2/s)  is  similar  to  that  obtained  in  our  simulation  with  the  KC04  LCMP  (400  cm2/s)  and  by 
KC04  (430  cm2/s).  The  shape  of  the  original  MYL2.5  and  KC04  LCMP  vertical  mixing  coefficient 
profiles  in  Fig.  2  is  more  similar  to  those  in  MW97  than  in  KC04  in  that  the  profiles  in  Fig.  2  and 
in  MW97  have  a  roughly  parabolic  shape,  whereas  the  profiles  in  KC04  reduce  fairly  abruptly  from 
a  maximum  value  to  a  value  close  to  zero  within  a  short  distance  near  the  base  of  the  SML. 

For  the  H15  LCMP,  the  maximum  value  of  Km  is  about  1500  cm2/s  at  7  m  depth.  This  is 
almost  4  times  larger  than  the  maximum  value  for  the  KC04  LCMP  and  about  7.5  times  larger 
than  the  maximum  value  for  the  original  MYL2.5  MLM.  This  large  value  of  Km  for  H15  is  due 
to  large  values  computed  for  H15’s  TKE  (see  Fig.  2)  and  stability  functions,  e.g.,  the  value  of  Sm 
corresponding  to  the  maximum  value  of  I\m  for  H15  is  about  1.15  compared  with  a  value  of  0.43 
for  KC04. 

The  component  of  the  model’s  Eulerian  velocity  aligned  with  the  wind  u  in  Fig.  2  shows  reduced 
shear  with  the  stronger  LC  mixing  and,  as  a  result,  the  value  of  u  at  the  surface  is  noticeably 
reduced,  which  is  also  shown  in  KC04  (their  Fig.  1)  and  MW97  (their  Fig.  2).  The  u  profile  for 
H15  in  Fig.  2  shows  less  shear  than  for  KC04  and  the  value  at  the  surface  is  reduced  as  well  to  the 
point  that  it  is  upwind.  Note  that  the  velocity  profiles  in  our  Fig.  2,  in  Fig.  2  in  MW97,  and  in 
Fig.  1  in  KC04  are  for  just  the  models’  Eulerian  current  and  do  not  include  the  SDC. 

With  the  inclusion  of  the  SDC  in  the  Coriolis  term  in  the  momentum  equations,  the  net 
transport  of  u  is  not  zero,  but  balances  the  net  transport  of  the  SDC  us,  so  that  the  combined  net 
downwind  transport  of  the  Eulerian  current  plus  the  SDC  is  zero.  Hence,  since  the  net  transport 
of  the  SDC  is  downwind  for  this  test  problem,  the  net  transport  of  the  downwind  component  of 
the  Eulerian  current  in  the  models  is  upwind  as  shown  in  Fig.  2  and  by  KC04  and  MW97. 

A  difference  in  MW97  from  the  simulation  here  and  in  KC04  for  the  case  without  the  LC 
mixing  is  that  MW97  drop  the  SDC  from  the  simulation  completely,  whereas  here  and  in  KC04 
the  SDC  is  retained  in  the  calculation  of  the  Coriolis  term  and  is  only  turned  off  for  the  calculation 
of  the  vertical  mixing  in  the  MYL2.5  turbulence  model.  Hence,  in  Fig.  2  and  in  KC04’s  Fig.  1, 
the  net  transport  of  u  for  the  case  without  LC  mixing  is  still  upwind  to  balance  the  net  downwind 
transport  of  the  SDC,  whereas  in  MW97,  the  net  transport  of  u  for  the  case  without  LC  mixing  is 
zero  (their  Fig.  3). 

The  cross-wind  velocity  v  in  Fig.  2  also  shows  reduced  shear  with  the  LC  mixing,  similar  to 
the  cross-wind  velocity  in  MW97  (their  Fig.  2),  and  the  surface  value  of  the  cross-wind  velocity  is, 
as  a  result,  reduced.  The  reduction  of  the  vertical  shear  and  the  surface  value  of  the  cross-wind 
velocity  is  greater  for  H15  than  for  KC04  because  of  H15’s  stronger  mixing.  The  net  transport  of 
the  cross-wind  velocity  must  equal  the  Ekman  transport  for  all  of  these  simulations  once  they  have 
reached  near-equilibrium. 

The  profiles  in  Fig.  2  indicate  slightly  deeper  mixing  for  the  cases  with  the  LCMP,  i.e.,  for  the 
case  without  the  LCMP,  the  MLD  increases  from  the  initial  value  of  33  to  34  m,  and  for  the  cases 
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with  the  LCMP  the  MLD  increases  from  33  to  about  37  m  for  both  KC04  and  H15.  A  deeper  MLD 
with  the  LC  mixing  is  also  indicated  in  the  results  shown  by  MW97  and  KC04. 

The  velocity  profiles  in  KC04  (their  Fig.  1)  do  not  appear  completely  consistent,  since  the 
transports  for  the  cases  with  and  without  LC  mixing  should  be  the  same  and  they  are  somewhat 
different.  This  may  be  because  the  velocity  profiles  shown  by  KC04  have  not  reached  equilibrium 
and  there  is  some  residual  inertial  motion  present  in  the  velocity  profiles  that  they  show. 

In  summary,  the  results  obtained  for  the  KC04  LCMP  for  this  simple,  wind-mixing  test  case 
are  generally  similar  to  the  results  reported  in  KC04  and  MW97.  The  values  of  the  maximum 
mixing  rates  for  momentum  are  roughly  doubled  by  including  the  LCMP  in  the  results  here  and  in 
KC04.  This  is  less  than  the  factor  of  three  increase  obtained  by  MW97,  but  the  maximum  mixing 
rate  for  momentum  obtained  here  with  KC04’s  LCMP  of  400  cnr/s  is  similar  to  that  obtained  by 
KC04  and  MW97  with  the  LC  mixing.  The  approximately  parabolic  shape  of  the  eddy  coefficient 
mixing  profiles  obtained  here  with  NCOM  are  roughly  similar  to  the  shape  of  the  eddy  coefficient 
profiles  obtained  by  MW97,  but  are  somewhat  different  from  the  shape  of  the  profiles  in  KC04. 

For  the  H15  LCMP,  the  TKE  is  higher  than  for  the  KC04  LCMP,  and  the  TKE  decreases  in 
the  upper  3  m,  whereas  the  TKE  for  the  KC04  LCMP  and  the  original  MYL2.5  MLM  increase 
near  the  surface  due  to  the  surface  flux  of  TKE  from  the  waves.  The  vertical  mixing  for  H15  is 
much  stronger  than  for  KC04,  i.e.,  the  maximum  value  of  Km  for  H15  is  1500  cm2/s  vs  400  cm2/s 
for  KC04.  The  stronger  mixing  for  the  H15  LCMP  results  in  less  vertical  shear  in  the  wind-driven 
surface  velocity  profiles  than  for  the  KC04  LCMP. 

4.  TEST  OF  LC  MIXING  IN  MYL2.5  MLM  AT  OWS  PAPA 

4.1  Description  of  Papa  Simulations 

OWS  Papa  was  located  in  the  northeast  Pacific  at  about  145°W,  50°N  from  1949  to  1981.  The 
availability  of  long  time  series  of  meteorological  and  ocean  subsurface  measurements  at  Papa,  and 
the  fact  that  the  effects  of  advection  on  the  heat  budget  tend  to  be  small  in  this  area,  have  long 
made  Papa  a  popular  location  to  test  and  evaluate  upper-ocean  MLMs  (Denman  and  Miyake  1973; 
Mellor  and  Durbin  1975;  Martin  1985;  Martin  1986;  Gaspar  1988;  Large  et  al.  1994;  Kantha  and 
Clayson  1994;  Large  1996). 

Simulations  were  conducted  at  OWS  Papa  for  the  year  1961  using  NCOM  run  in  pseudo  one¬ 
dimensional  mode  and  using  the  MYL2.5  turbulence  model  both  with  and  without  the  LCMPs  of 
KC04  and  H15. 

For  the  simulations  at  OWS  Papa,  a  vertical  grid  of  100  layers  was  used,  with  a  surface-layer 
thickness  of  1  nr  and  a  smooth  stretching  to  a  maximum  depth  of  5500  m.  Hence,  each  layer  is 
thicker  than  the  layer  above  by  about  6%.  This  is  higher  vertical  resolution  than  is  typically  used 
with  NCOM.  However,  sensitivity  tests  have  shown  that  the  Mellor- Yamada-type  mixing  schemes 
in  NCOM  tend  to  predict  a  slighly  deeper  mixed  layer  and  a  smoother,  better-resolved  deepening 
of  the  mixed  layer,  as  the  vertical  resolution  is  increased;  this  is  especially  noticeable  at  OWS  Papa 
when  the  mixed  layer  deepens  in  the  fall  (Martin  and  Hogan,  2013).  Hence,  high  vertical  grid 
resolution  was  used  for  the  Papa  simulations  to  better  illustrate  the  performance  capabilities  of  the 
MLMs  in  comparisons  with  the  observations  at  Papa. 
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Hourly  surface  wind  stresses  and  heat  fluxes  were  computed  from  the  3-hourly  surface  marine 
observations  at  Papa  using  fairly  standard  formulas  as  discussed  in  Martin  (1985,  1986).  Note 
that  all  the  heat  fluxes  were  pre-computed  using  the  surface  marine  observations  from  Papa,  i.e. , 
the  NCOM-predicted  SST  was  not  used  in  the  calculation  of  any  of  the  heat  fluxes.  The  mean 
computed  heat  fluxes  over  the  year  were  105  W/m2  for  the  solar  radiation,  -88  W/m2  for  the  surface 
heat  loss  due  to  the  net  longwave  and  latent  and  sensible  heat  fluxes,  and  18  W/m2  for  the  net 
heat  flux.  These  pre-computed  heat  fluxes  are  fairly  consistent  with  the  seasonal  changes  in  the 
observed  upper-ocean  heat  content  at  Papa  during  the  spring  and  summer  of  1961  (Martin  et  al. 
2012),  and  this  is  the  reason  that  all  the  heat  fluxes  are  fixed  and  none  (e.g.,  the  latent  and  sensible 
heat  fluxes)  are  computed  using  the  model-predicted  SST. 

The  solar  extinction  was  taken  to  be  Jerlov  Type  II  (Jerlov  1968).  The  solar  extinction  (i.e.,  the 
depth  of  penetration  of  the  solar  radiation)  can  significantly  affect  the  developement  of  the  upper- 
ocean  thermal  structure  in  the  summer  (Martin,  1985),  however,  this  effect  tends  to  be  larger  at 
lower  latitudes  where  the  difference  in  the  summer /winter  heating  is  less  pronounced  and  the  winds 
and  the  deepening  of  the  SML  in  the  summer  by  the  winds  are  reduced  (Martin  and  Allard,  1993). 

The  ambient /background  viscosity  and  diffusivity  below  the  SML  were  taken  to  be  0.02  cm2/s. 
The  ambient  diffusivity  affects  the  predicted  SST  in  the  summer  at  Papa,  since  a  larger  value  of 
the  diffusivity  causes  a  larger  downward  transport  of  heat  out  of  the  mixed  layer  and  into  the 
thermocline  during  the  summer.  The  sharp  seasonal  thermocline  that  develops  at  Papa  during  the 
heating  season  suggests  that  the  ambient  diffusivity  below  the  SML  at  this  time  is  fairly  small. 
The  relative  magnitudes  of  the  ambient  diffusivity  and  viscosity  affect  the  deepening  of  the  mixed 
layer,  e.g.,  increasing  the  ambient  viscosity  relative  to  the  ambient  diffusivity  tends  to  reduce  the 
vertical  shear  of  the  horizontal  velocity  near  the  base  of  the  SML  and  increase  the  local  Richardson 
number  and,  hence,  reduce  the  turbulent  mixing  (and  visa  versa). 

A  linear  damping  term  in  the  momentum  equations  with  a  damping  time  scale  of  10  d  was 
used  to  provide  some  damping  of  the  inertial  oscillations,  since  the  normal  horizontal  and  vertical 
dissipation  of  these  motions  cannot  be  accounted  for  in  one-dimensional  simulations  such  as  those 
being  conducted  here.  Without  such  damping,  inertial  oscillations  generated  in  the  winter  and 
spring  when  the  SML  is  relatively  deep  can  persist  below  the  SML  in  the  summer  when  the  SML 
is  shallow,  and  can  affect  the  deepening  of  the  SML  in  the  fall  when  these  persisting  inertial 
oscillations  are  encountered  by  the  deepening  SML.  Although  interactions  between  a  deepening 
SML  and  existing  inertial  oscillations  commonly  occur  in  the  real  ocean,  inertial  oscillations  in  the 
ocean  are  governed  by  three-dimensional  processes  that  are  not  accounted  for  in  a  one-dimensional 
simulation,  and  it  was  considered  best  to  avoid  the  artificial  occurrance  of  encountering  inertial 
oscillations  generated  and  (unrealistically)  left  behind  in  the  deeper  water  weeks  or  months  earlier. 

The  initial  condition  for  the  simulations  at  Papa  was  a  temperature  profile  from  the  bathyther¬ 
mograph  (BT)  observations  at  Papa  for  1  January  1961,  and  a  climatological  salinity  profile  for  the 
location  of  Papa  for  the  beginning  of  January  from  the  Levitus  (1982)  climatology.  The  salinity 
profiles  at  Papa  feature  a  strong  halocline  between  about  100  and  200  nr  that  helps  limit  the  depth 
of  mixing  in  the  winter. 

Including  the  effect  of  the  LCMPs  of  KC04  and  H15  on  the  simulation  of  the  SML  at  Papa 
requires  an  estimate  of  the  SDC.  Since  we  don’t  have  wave  observations  for  Papa  for  1961,  the 
SDC  profile  was  estimated  from  the  observed  surface  winds  at  Papa.  As  discussed  in  Martin  et  al. 
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(2013),  there  are  a  number  of  ways  this  can  be  done,  e.g.,  (1)  compute  the  SDC  from  the  wind 
stress  using  empirical  functions  like  those  used  by  KC04  for  their  simulations  at  OWS  Papa,  (2) 
compute  the  SDC  from  the  wind  speed  using  wave  energy  spectra  time-averaged  for  different  wind 
speeds  using  wind  and  wave  data  from  NOAA  Buoy  46005  located  in  the  northeast  Pacific  about 
1100  km  east-southeast  of  Papa,  and  (3)  compute  the  SDC  from  the  waves  predicted  by  a  wave 
model  forced  by  the  observed  winds  at  Papa  for  1961.  Martin  et  al.  (2013)  found  that  these  three 
methods  for  estimating  waves  at  Papa  resulted  in,  from  (1)  to  (3),  successively  smaller  wave  effects. 
The  method  used  in  this  report  is  (2),  i.e.,  to  estimate  the  waves  from  the  wind  based  on  relations 
between  the  wind  and  waves  derived  from  observations  from  NOAA  Buoy  46005.  One  reason  for 
this  choice  is  that  this  method  includes  swell,  which  always  seems  to  be  present  in  this  area,  even 
when  the  winds  are  light  (Martin  et  al.  2012). 

Nine  years  of  wind  and  wave  data  from  NOAA  Buoy  46005  in  the  northeast  Pacific  were  used 
to  compute  average  wave  spectra  within  1-m/s-wide,  wind-speed  bins.  These  average  wave  spectra 
were  then  used  to  compute  hourly  SDC  profiles  at  OWS  Papa  for  1961  based  on  the  hourly  wind 
speed  at  Papa.  The  details  of  this  calculation  can  be  found  in  Appendix  B.  Buoy  46005  is  located 
at  131.001°W,  46.100°N,  which  is  about  1100  km  ESE  of  the  location  of  OWS  Papa.  Hence,  the 
wind  and  wave  conditions  in  the  area  of  the  buoy  should  be  roughly  similar  to  those  near  OWS 
Papa.  This  calculation  of  the  SDC  at  Papa  will  be  referred  to  as  the  Buoy  Wave  Spectra  (BWS) 
SDC. 

4.2  Papa  Simulations  with  MYL2.5  MLM  without  LC  mixing 

Figure  3  shows  results  from  a  simulation  of  the  SML  at  Papa  for  1961  using  just  the  MYL2.5 
MLM  with  no  LCMP.  However,  note  that  the  SDC  is  used  in  the  Coriolis  terms  of  the  ocean  model’s 
momentum  equations  as  discussed  in  Section  3.  In  the  figure,  the  results  of  the  simulation  (red)  are 
compared  with  the  observations  at  Papa  (black).  The  Helds  compared  are  SST,  MLD,  and  isotherm 
depths  (ISODs).  The  MLD  shown  in  the  figure  is  the  depth  at  which  the  temperature  becomes 
0.2°C  less  than  the  SST.  The  isotherm  depths  are  plotted  at  1°C  intervals  and  are  temporally 
filtered  to  reduce  the  high-frequency  variability  and  make  them  easier  to  discerne. 

As  in  Martin  (1985,  1986),  the  SST  is  too  high  in  summer  and,  consistent  with  this,  the  MLD 
and  ISODs  are  too  shallow.  In  Fig.  3,  the  summer  SST  is  about  2.5°C  too  high  and  the  MLD  and 
ISODs  are  10  to  20  m  too  shallow  in  August  and  September.  The  agreement  in  Fig.  3  is  a  bit  worse 
than  in  Martin  (1985,  1986)  because,  for  this  simulation,  the  net  surface  heat  flux  has  been  set 
slightly  higher  to  provide  more  net  heating  and  better  agreement  with  the  observed  upper-ocean 
heat  content  at  Papa  in  the  summer  (Martin  et  al.  2012). 

There  are  various  methods  of  enhancing  the  depth  of  mixing  obtained  with  the  MYL2.5  MLM 
(Martin  et  al.  2013).  KC94  used  a  parameterization  of  shear  instability  computed  from  the  local 
Richardson  number  that  was  taken  from  Large  et  al.  (1994,  hereafter  referred  to  as  L94).  The 
L94  shear-instability  mixing  parameterization  (SIMP)  provides  a  moderate  level  of  mixing  for 
Richardson  numbers  above  the  normal  turbulence  cutoff  of  0.2-0.25  up  to  a  Richardson  number  of 
0.7.  The  maximum  mixing  rate  for  the  L94  SIMP  is  50  cm2/s  at  a  Richarson  number  of  zero  that 
decreases  to  zero  at  a  Richardson  number  of  0.7.  A  more  complete  description  of  the  L94  SIMP  is 
provided  in  Appendix  C. 

KC94  described  the  parameterization  as  accounting  for  intermittent  mixing  due  to  reduction  of 
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Fig.  3  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  MYL2.5  MLM  with 

no  LCMP  (red).  ISODs  are  temporally  filtered. 
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the  local  Ri  by  internal  waves  below  the  SML.  L94  described  this  parameterization  as  accounting 
for  shear  instability  due  to  resolved  vertical  shear  below  the  SML  (and  they  provided  a  separate 
parameterization  for  mixing  by  internal  waves).  Hence,  L94  and  KC94  do  not  seem  to  be  in 
complete  agreement  as  to  what  the  L94  SIMP  represents;  however,  its  main  effect  is  to  provide 
some  additional  mixing  for  Richardson  numbers  above  the  normal  critical  value,  which  effectively 
increases  the  predicted  depth  of  the  SML.  Both  L94  and  KC94  reported  that  the  L94  SIMP  improved 
agreement  with  observed  MLDs  for  several  open  ocean  data  sets,  including  data  from  OWS  Papa. 
Hence,  the  L94  SIMP  may  represent  some  unresolved  or  unaccounted  for  processes,  including 
perhaps  LC,  that  contribute  to  the  deepening  of  the  SML.  The  L94  SIMP  has  been  available  to 
enhance  the  vertical  mixing  of  the  Mellor-Yamada  Level  2  turbulence  model  in  NCOM,  and  NCOM 
was  recently  modified  to  allow  the  L94  SIMP  to  be  used  with  the  MYL2.5  turbulence  model  as 
well. 

Figure  4  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  using 
the  MYL2.5  MLM  with  the  L94  SIMP.  With  the  L94  SIMP,  the  depth  of  the  SML  is  significantly 
increased  in  the  summer  and  fall  and  the  MLD,  ISODs,  and  SST  all  show  better  agreement  with 
the  observed  values.  That  the  SST  and  the  MLD  and  ISODs  simultaneously  show  fairly  good 
agreement  with  the  observations  suggests  that  the  net  surface  heat  flux  used  for  the  simulation  is 
a  fairly  good  match  for  the  actual  changes  in  the  upper-ocean  heat  content  at  Papa  for  1961. 

4.3  Papa  Simulations  with  KC04  parameterization  of  LC  mixing 

Figure  5  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  using 
the  MYL2.5  MLM  with  the  KC04  LCMP.  Figure  5  shows  that  the  simulation  is  slightly  improved 
relative  to  that  for  the  MYL2.5  MLM  without  the  LCMP  in  Fig.  3,  but  the  SST  is  still  about  1.5°C 
too  high  in  August  and  September  and,  consistent  with  this,  the  MLD  and  ISODs  are  too  shallow. 
Hence,  the  use  the  KC04  LCMP  slightly  increases  the  predicted  MLD,  but  the  increase  does  not 
seem  to  be  sufficient  to  account  for  the  observed  MLD  at  Papa  in  the  summer. 

Figure  6  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  using 
the  MYL2.5  MLM  with  both  the  KC04  LCMP  and  the  L94  SIMP.  As  for  the  simulations  with 
the  MYL2.5  MLM  in  Fig.  3  and  4,  the  addition  of  the  L94  SIMP  significantly  increases  the  MLD 
during  the  summer  and  fall.  The  MLD  and  ISODs  at  this  time  are  now  slightly  too  deep  and, 
consistent  with  this,  the  SST  is  slightly  too  cool.  Note  that  this  simulation  could  be  adjusted  to 
better  agree  with  the  observations  by  reducing  the  value  of  the  Richardson  number  at  which  the 
L94  SIMP  cuts  off  mixing. 

Figure  7  shows  hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients 
Km  hi  the  SML  for  the  Papa  simulations  with  and  without  the  KC04  LCMP  that  are  shown  in 
Fig.  6  and  4,  respectively  (both  these  simulations  use  the  L94  SIMP).  Plots  of  the  ratio  of  both 
unfiltered  and  temporally  filtered  Km  values  are  shown  in  Fig.  7. 

The  unfiltered  plot  shows  that  the  range  of  the  ratio  of  the  Km  values  is  usually  between 
about  0.5  and  3.5,  with  occasional  lower  and  higher  values.  Values  less  than  one  indicate  that 
the  maximum  value  of  Km  for  the  simulation  with  the  KC04  LCMP  is  lower  than  that  for  the 
simulation  without  the  LCMP. 

The  plot  of  the  ratio  of  the  temporally  filtered  I\m  values  shows  a  range  between  about  1  and  3 
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Fig.  4  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  MYL2.5  MLM  with 

the  L94  SIMP  (red).  ISODs  are  temporally  filtered. 
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Fig.  5  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  KC04  LCMP  (red). 

ISODs  are  temporally  filtered. 
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Fig.  6  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  using  both  the  KC04  LCMP 

and  the  L94  SIMP  (red).  ISODs  are  temporally  filtered. 
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Fig.  7  —  Hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients  Km  in  the  SML  for 
Papa  simulations  with  and  without  the  KC04  LCMP  (both  include  the  L94  SIMP).  The  left  plot  shows  the  ratio  of 
unfiltered  Km  values  and  the  right  plot  shows  the  ratio  of  temporally  filtered  Km  values. 
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Fig.  8  —  Surface  current  magnitude  for  simulation  without  the  KC04  LCMP  (left)  and  with  the  KC04  LCMP  (right), 

both  plots  have  had  a  strong  temporal  filter  applied. 


with  occasional  higher  values.  The  ratio  is  smaller  in  the  fall  and  winter  when  free  convection,  due 
mainly  to  surface  cooling,  plays  a  larger  role  in  the  vertical  mixing,  and  is  larger  in  the  spring  and 
summer  when  forced  convection  due  to  wind  forcing  and  the  resulting  shear  production  of  TKE 
plays  a  larger  role.  The  mean  value  of  the  ratio  in  the  spring  and  summer  is  about  2.3. 

Figure  8  shows  hourly  values  of  the  magnitude  of  the  surface  current  for  the  Papa  simulations 
without  and  with  the  KC04  LCMP  (note  that  these  time  series  have  been  temporally  filtered). 
The  mean  value  of  the  magnitude  of  the  surface  current  is  reduced  about  10%  over  the  year  when 
the  KC04  LCMP  is  used.  This  is  somewhat  less  than  the  reduction  of  the  surface  current  by 
the  KC04  LCMP  shown  in  Fig.  2.  This  is  partly  because  of  the  frequent  occurrence  of  inertial 
oscillations  in  the  Papa  simulations,  which  tend  to  be  fairly  uniform  with  depth  in  the  SML  and 
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Fig.  9  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  H15  LCMP  (red). 

ISODs  are  temporally  filtered. 


are  not  significantly  affected  by  the  strength  of  the  vertical  mixing,  whereas  for  the  velocity  profiles 
in  Fig.  2  the  inertial  oscillations  were  suppressed  by  turning  on  the  wind  forcing  gradually  using  a 
linear  ramp  over  one  inertial  period. 

4.4  Papa  Simulations  with  H15  parameterization  of  LC  mixing 

Figure  9  shows  results  from  a  simulation  of  the  upper-ocean  thermal  structure  at  Papa  using 
the  MYL2.5  MLM  with  the  H15  LCMP.  Note  that  this  simulation  is  not  using  the  L94  SIMP. 
Comparison  of  this  simulation  with  the  earlier  Papa  simulations  indicates  that  the  H15  LCMP 
predicts  an  overall  deeper  SML  than  the  original  MYL2.5  MLM  and  the  KC04  LCMP  (Fig.  3  and 
5,  respectively),  and  provides  a  good  simulation  of  the  SST,  MLD,  and  ISODs  at  Papa  without 
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much  need  for  additional  mixing  enhancements.  The  SST  is  slightly  too  warm  in  summer  and  the 
MLD  and  IDODs  are,  correspondingly,  slightly  too  shallow. 

Since  the  H15  LCMP  has  shear-production  terms  in  the  TKE  and  TLS  equations  (19-20)  that 
depend  only  on  the  shear  produced  by  the  SDC,  there  was  some  concern  whether  the  H15  LCMP 
would  predict  adequate  shallowing  during  light- wind  conditions  at  Papa,  since  some  swell  is  always 
present.  However,  Fig.  9  indicates  that  the  H15  LCMP  predicts  shallowing  and  SST  peaks  during 
calm  wind  conditions  similar  to  the  original  MYL2.5  MLM  and  the  KC04  LCMP. 

Figure  10  shows  hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coeffi¬ 
cients  Km  in  the  SML  for  the  H15  and  KC04  Papa  simulations  shown  in  Fig.  9  and  6,  respectively 
(the  KC04  simulation  here  uses  the  L94  SIMP  but  the  H15  simulation  does  not).  Plots  of  the  ratio 
of  both  unfiltered  and  temporally  filtered  Km  values  are  shown.  The  unfiltered  Km  ratios  show 
a  large  range.  The  filtered  Km  ratios  show  values  between  about  1  and  5,  with  a  mean  value  in 
the  spring  and  summer  of  2.8.  Hence,  as  in  the  simple  mixing  test  case  (Fig.  2),  the  H15  LCMP 
predicts  significantly  higher  mixing  rates  for  the  shear-generated  TKE  than  the  KC04  LCMP.  A 
plot  of  the  time-filtered  surface  current  magnitudes  for  these  two  simulations  (not  shown)  indicates 
that  the  mean  surface  currents  are  similar. 

Li  et  al.  (2005)  describe  a  regime  diagram  for  classifying  turbulent  mixing  in  the  upper  ocean. 
Such  a  diagram  is  presented  here  in  Fig.  11  for  OWS  Papa  for  1961  using  the  SDC  estimated  from 
the  wind  and  wave  data  from  NOAA  Buoy  46005  (i.e.,  the  BWS).  The  figure  shows  a  scatter  plot 
of  hourly  values  of  the  Hoenikker  number  versus  the  turbulent  Langmuir  number  for  Papa  for  1961. 
The  Hoenikker  number  H0  is  given  by 


H0 


4Bods 
14(0 )rt*  ’ 


(37) 


where  Bq  is  the  total  surface  buoyancy  flux,  ds  is  the  e-folding  depth  of  the  SDC,  14(0)  is  the  surface 
value  of  the  SDC,  and  u*  is  the  water-side  surface  friction  velocity.  The  turbulent  Langmuir  number 
Lat  is  given  by 


Lat 


f  u*  \  1/2 
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(38) 


The  diagram  is  divided  into  three  regimes  according  to  the  type  of  turbulent  mixing:  Langmuir, 
convection,  or  shear,  that  dominates  within  that  region.  The  boundaries  of  the  three  regimes 
were  determined  by  Li  et  al.  (2005)  based  on  LES  experiments.  The  regimes  are  characterized  by 
the  relative  magnitudes  of  the  components  of  the  turbulent  mixing,  with  the  vertical  component 
dominating  in  convective  mixing,  the  downwind  component  dominating  in  shear  mixing,  and  the 
cross-wind  and  vertical  components  dominating  the  downwind  component  in  Langmuir  mixing. 


Figure  11  indicates  that  the  upper-ocean  mixing  at  Papa  for  1961  is  generally  in  the  Langmuir 
and  convection  regimes.  Only  a  few  points  occur  in  the  shear  regime,  and  these  occur  close  to  the 
boundary  between  the  shear  and  Langmuir  regimes.  This  is  consistent  with  the  finding  of  Li  et  al. 
(2005)  that  wind  mixing  in  the  upper-ocean  typically  occurs  in  the  Langmuir  regime. 

As  noted  in  Section  2.3,  the  behavior  of  the  H15  turbulence  parameterization  when  the  SDC 
is  zero  is  of  interest,  since  the  SDC  in  the  SML  will  be  small  when  the  waves  are  small,  and  will  be 
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Fig.  10  —  Hourly  values  of  the  ratio  of  the  maximum  value  of  the  vertical  mixing  coefficients  Km  in  the  SML  for 
Papa  simulations  with  the  H15  LCMP  and  with  the  KC04  LCMP.  The  left  plot  shows  the  ratio  of  unfiltered  Km 
values  and  the  right  plot  shows  the  ratio  of  temporally  filtered  Km  values. 
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Hoen i kker  Number  vs  Turbulent  Langmuir  Number 


LoglO  (Turbulent  Langmuir  Number) 

Fig.  11  —  Turbulent  regime  diagram  for  OWS  Papa  for  1961  using  the  SDC  computed  from  the  BWS.  The  points 
plotted  are  hourly  values  of  the  loglO  of  the  Hoenikker  number  (with  an  offset  of  0.01)  versus  the  loglO  of  the  turbulent 
Langmuir  number. 
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Fig.  12  —  SST,  MLD,  and  ISODs  at  OWS  Papa  for  1961:  observed  (black),  simulated  with  the  H15  LCMP  with 

SDC=0  (red).  ISODs  are  temporally  filtered. 

zero  in  deep  mixed  layers  below  the  influence  of  the  surface  waves.  Figure  12  shows  results  from 
a  simulation  at  Papa  with  H15  with  the  SDC  set  equal  to  zero.  The  MLD  is  shallower  than  in 
Fig.  9  for  H15  with  the  SDC,  and  is  similar  to  the  MLD  for  the  MYL2.5  turbulence  model  without 
LCMP  mixing  in  Fig.  3. 

5.  TEST  OF  LC  MIXING  IN  MYL2.5  MLM  FOR  HURRICANE  IVAN  SIMULA¬ 
TIONS 

Simulations  were  run  with  NCOM  in  the  Gulf  of  Mexico  (GoM)  with  wave  forcing  from  SWAN 
to  look  at  the  effect  of  the  KC04  and  H15  LCMPs  on  a  simulation  of  the  response  of  the  ocean 
in  the  GoM  to  Hurricane  Ivan  in  September  2004.  Atmospheric  forcing  was  provided  from  a 
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simulation  of  Hurricane  Ivan  with  the  COAMPS  atmospheric  model.  COAMPS,  which  stands 
for  Coupled  Ocean/ Atmosphere  Mesoscale  Prediction  System  (Chen  et  al.  2010),  is  currently  a 
coupled  atmosphere-ocean-wave  modeling  system.  However,  for  the  simulations  conducted  here, 
the  atmosphere,  ocean,  and  wave  models  were  run  separately.  Within  the  rest  of  this  report, 
COAMPS  will  be  used  to  refer  to  just  the  COAMPS  atmospheric  model  (Hodur  1997). 

NCOM  was  run  with  a  horizontal  grid  resolution  of  about  4  km.  NCOM’s  vertical  grid  con¬ 
sisted  of  40  layers  with  a  surface  layer  thickness  of  1  m  and  a  smooth  stretching  to  a  maximum 
depth  of  5500  m.  Initial  and  boundary  conditions  for  NCOM  were  provided  by  the  Global  NCOM 
model  that  is  run  operationally  at  the  Naval  Oceanographic  Office  at  Stennis  Space  Center,  Mis¬ 
sissippi  (Barron  et  al.  2004).  The  COAMPS  atmospheric  model  was  run  as  a  triply- nested  system 
with  grid  resolutions  of  81,  27,  and  9  km.  Initial  and  boundary  conditions  for  COAMPS  were 
provided  by  the  Navy  Operational  Global  Atmospheric  Prediction  System  (NOGAPS),  which  was 
run  at  the  Navy’s  Fleet  Numerical  Meteorology  and  Oceanography  Center  (FNMOC)  in  Monterey, 
California  (Rosmond  et  al.  2002).  Hourly  surface  fields  from  COAMPS’s  innermost  grid  were  used 
to  force  SWAN  and  NCOM.  These  fields  include  the  surface  winds  for  SWAN,  and  the  surface  wind 
stress,  surface  pressure,  solar  radiation,  net  longwave  radiation,  latent  and  sensible  heat  fluxes, 
evaporation,  and  precipitation  for  NCOM.  SWAN  was  run  on  the  same  4-krn  horizontal  grid  as 
NCOM. 

Wave  fields  for  forcing  NCOM  were  output  hourly  from  SWAN  and  include  the  surface  wave 
radiation  stress,  the  SDC  profiles,  and  the  characteristic  amplitude,  frequency,  and  direction  of 
the  orbital  wave  motion  near  the  ocean  bottom.  The  wave  radiation  stress  forcing  in  NCOM 
is  implemented  similar  to  the  surface  wind  stress  as  a  surface  stress.  In  the  open  sea,  the  wave 
radiation  stress  forcing  of  the  ocean  is  relatively  small,  i.e.,  generally  less  than  10%  of  the  magnitude 
of  the  wind  stress.  The  SDC,  as  noted  in  Section  2,  is  used  in  NCOM’s  Coriolis  and  advection  terms, 
in  the  continuity  equation,  and  also  in  the  KC04  and  H15  LCMPs.  The  wave  orbital  motions  near 
the  bottom  act  to  enhance  the  ocean  model’s  bottom  drag  in  areas  where  the  bottom  is  sufficiently 
shallow  to  feel  the  orbital  wave  motions  (Grant  and  Masden  1979;  Soulsby  1995).  However,  the 
focus  here  is  only  on  the  effect  of  the  KC04  and  H15  LCMPs;  hence,  the  NCOM  simulations  that  are 
compared  all  have  full  wave  forcing  from  SWAN,  and  the  only  difference  between  them  is  whether 
the  LCMPs  of  KC04  or  H15  are  used  or  not. 

The  ocean  model  simulation  was  started  at  00Z  on  12  September  2004,  at  which  time  Hurricane 
Ivan  was  in  the  central  Caribbean  Sea.  Figure  13  shows  the  SST  and  sea-surface  velocity  (SSV) 
vectors  for  the  simulation  with  the  H15  LCMP  for  09Z  15  September.  Hurricane  Ivan  is  near  the 
center  of  the  GoM  at  this  time.  The  surface  wind  vectors  from  the  COAMPS  atmospheric  model 
are  also  shown  in  the  plot.  Note  that  the  surface  current  vectors  in  Fig.  13  include  both  the  NCOM 
surface  current  and  the  SDC  from  SWAN.  The  hurricane  winds  generate  strong  cooling  along  the 
hurricane  track,  with  the  strongest  cooling  usually  occurring  just  to  the  right  of  the  track  (Price 
1981).  The  equivalent  plot  for  the  simulation  with  the  KC04  LCMP  looks  very  similar  (this  plot  is 
shown  in  M13  in  Fig.  14,  but  is  not  shown  here). 

Figure  14  shows  the  difference  in  the  SST  between  the  simulation  with  the  H15  LCMP  (Fig.  13) 
and  with  no  LCMP  (not  shown)  at  09Z  15  September  2004.  The  simulation  with  the  H15  LCMP 
generally  shows  more  cooling  of  the  SST  near  the  hurricane,  as  indicated  by  the  areas  of  negative 
temperature  difference;  the  difference  is  typically  a  few  tenths  of  a  °C  and  the  maximum  difference 
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Fig.  13  —  SST  (color  contours),  SSV  (black  vectors),  and  surface  winds  (white  vectors)  at  09Z  15  September  2004 

for  simulation  of  Hurricane  Ivan  with  the  H15  LCMP 
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Fig.  14  —  Difference  in  SST  for  Hurricane  Ivan  simulations  with  the  H15  LCMP  and  without  a  LCMP  at  09Z  15 

September  2004. 


is  about  0.8°C.  The  lower  SST  for  the  simulation  with  the  H15  LCMP  indicates  deeper  mixing, 
which  is  consistent  with  the  results  at  OWS  Papa  in  Section  4. 

Figure  15  shows  the  equivalent  plot  to  Fig.  14  for  the  KC04  LCMP  (note,  however,  that  the 
temperature  range  in  Fig.  15  is  half  that  used  for  Fig.  14).  The  additional  cooling  for  the  KC04 
LCMP  compared  to  the  original  MYL2.5  MLM  without  an  LCMP  is  generally  small,  i.e.,  less  than 
0.5  °C,  and  is  significantly  less  than  that  in  Fig.  14,  which  is  consistent  with  the  results  at  OWS 
Papa. 

Figure  16  shows  the  difference  in  the  SSV  between  the  simulation  with  the  H15  LCMP  and 
without  at  09Z  15  September  2004.  The  color  contours  show  the  difference  in  the  magnitude  of 
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Fig.  15  —  Difference  in  SST  for  Hurricane  Ivan  simulations  with  KC04  LCMP  and  with  no  LCMP  at  09Z  15  September 

2004. 
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Fig.  16  -  Difference  in  SSV  for  Hurricane  Ivan  simulations  with  the  H15  LCMP  and  for  the  original  MYL2.5  MLM 
without  a  LCMP  at  09Z  15  September  2004.  The  color  contours  show  the  difference  in  the  magnitude  of  the  surface 
current  and  the  black  vectors  show  the  vector  velocity  difference.  The  white  vectors  show  the  surface  winds. 

the  surface  current  and  the  black  vectors  show  the  vector  velocity  difference.  The  fact  that  the 
direction  of  the  velocity  difference  vectors  is  opposite  to  the  wind  direction  indicates  that  the 
stronger  mixing  with  the  H15  LCMP  reduces  the  wind-driven  surface  currents.  Near  the  hurricane, 
the  surface  currents  are  reduced  by  30-50  cm/s  or  more.  Some  areas  show  the  surface  current 
magnitude  with  the  H15  LCMP  to  be  larger  (as  indicated  by  the  color  contours  in  Fig.  16)  and 
this  is  because  the  direction  of  the  prevailing  currents  in  these  areas  is  nearly  opposite  to  the  wind 
direction,  and  so  a  reduction  in  the  downwind,  wind-driven  current  increases  the  total  current,  e.g., 
on  the  west  side  of  the  hurricane  there  is  an  anti-cylonic  Loop  Current  eddy  (see  Fig.  13). 

The  reduction  of  the  surface  current  for  the  KC04  LCMP  (Fig.  17)  is  less  than  for  the  H15 
LCMP  (Fig.  16),  since  the  vertical  mixing  rates  are  not  as  strong. 
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Fig.  17  —  Difference  in  SSV  for  Hurricane  Ivan  simulations  with  the  KC04  LCMP  and  without  a  LCMP  at  09Z  15 
September  2004.  The  color  contours  show  the  difference  in  the  magnitude  of  the  surface  current  and  the  black  vectors 
show  the  vector  velocity  difference.  The  white  vectors  show  the  surface  winds. 
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Fig.  18  -  -  Maximum  value  of  Km  in  SML  for  Hurricane  Ivan  simulation  with  H15  LCMP  at  09Z  15  September  2004. 

The  white  vectors  show  the  surface  winds. 

Figure  18  shows  the  maximum  value  of  the  vertical  mixing  coefficient  Km  in  the  SML  for  the 
simulation  with  the  H15  LCMP  at  09Z  15  September  2004.  The  Km  values  in  Fig.  18  were  spatially 
filtered  with  one  pass  of  a  3x3  box  filter.  The  strongest  mixing  near  the  hurricane  is  occurring  on 
the  northeast  side  of  the  eye,  with  maximum  values  of  Km  of  4  m2/s.  A  circular  area  of  weak 
mixing  occurs  under  the  eye  of  the  hurricane  due  to  the  lighter  winds  there. 

Figure  19  shows  the  maximum  value  of  the  vertical  mixing  coefficient  Km  in  the  SML  for  the 
simulation  with  the  KC04  LCMP  at  09Z  15  September  2004.  The  Km  values  in  Fig.  19  were  also 
spatially  filtered  with  one  pass  of  a  3x3  box  filter.  The  strongest  mixing  near  the  hurricane  is 
occurring  on  the  north  and  east  sides  of  the  eye,  with  maximum  values  of  Km  of  about  1  m2/s. 

The  local  time  for  Fig.  19  is  about  three  o’clock  in  the  morning;  hence,  there  is  relatively  strong 
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Fig.  19  —  Maximum  value  of  Km  in  SML  for  Hurricane  Ivan  simulations  with  KC04  LCMP  at  09Z  15  September 

2004.  The  white  vectors  show  the  surface  winds. 
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night-time  convective  mixing  occurring  in  some  areas,  e.g.,  in  the  western  GoM  and  the  northwest 
Caribbean.  On  the  afternoon  of  September  15  (local  time),  most  of  this  mixing  is  significantly 
reduced  due  to  day-time  solar  heating  (not  shown).  However,  the  mixing  near  the  hurricane  is 
primarily  wind-driven  and  remains  strong  through  the  day  and  night.  Also,  the  stronger  mixing 
south  of  and  just  north  of  Cuba  is  partly  due  to  this  area  being  on  the  warm  side  of  the  Loop 
Current  where  the  SML  is  generally  deeper  than  on  the  cold  side. 

Figure  20  shows  the  ratio  of  the  spatially-filtered  values  of  the  maximum  value  of  Km  in  the 
SML  for  the  simulations  with  the  H15  LCMP  and  without  a  LCMP  at  09Z  15  September  2004. 
Near  the  hurricane,  the  ratios  range  from  about  1  to  10.  This  indicates  that  the  H15  LCMP  is 
significantly  increasing  the  mixing  rates  in  the  areas  of  strong  hurricane  winds.  The  areas  where 
the  ratio  is  below  one  are  generally  occurring  in  areas  of  night-time  convection. 

Figure  21  shows  the  ratio  of  the  spatially- filtered  values  of  the  maximum  value  of  Km  in  the 
SML  for  the  simulations  with  the  KC04  LCMP  and  without  a  LCMP  at  09Z  15  September  2004. 
Near  the  hurricane,  the  ratios  range  from  1  to  4  and  are  frequently  in  the  range  of  2  to  3.  This 
indicates  that  the  KC04  scheme  is  significantly  increasing  the  mixing  rates  in  the  areas  of  strong 
hurricane  winds,  though  not  nearly  as  much  as  the  H15  LCMP.  The  areas  where  the  ratio  is  below 
one  are  generally  occurring  in  areas  of  night-time  convection. 

6.  SUMMARY 

Mellor-Yamada-type  turbulence  models  have  long  been  criticised  for  not  mixing  the  ocean’s 
surface  layer  strongly  or  deeply  enough  in  the  open  ocean.  Recent  LES  simulations  of  Langmuir 
circulation  in  the  upper  ocean,  e.g.,  MW97,  found  that  the  occurrence  of  LCs  significantly  increases 
the  net  rate  of  mixing  within  the  SML  and  slightly  increases  the  MLD.  The  original  Mellor-Yamada- 
type  turbulence  models  did  not  account  for  these  effects. 

Based  on  the  LES  results  of  MW97,  KC04  parameterized  and  tested  the  effects  of  enhanced 
mixing  by  LC  in  the  MYL2.5  turbulence  model.  KC04  parameterized  the  effect  of  LC  mixing  in 
the  MYL2.5  turbulence  model  by  adding  additional  shear-production  terms  to  the  TKE  and  TLS 
equations  of  the  turbulence  model.  These  additional  terms  consist  of  the  product  of  the  vertical 
turbulent  momentum  flux  and  the  vertical  shear  of  the  SDC  from  the  surface  wave  field. 

H15  noted  that  KC04  were  not  fully  consistent  in  their  application  of  the  CL  vortex  force 
to  the  ARSM  used  derive  the  MYL2.5  turbulence  model,  and  derived  alternative  modifications 
to  the  turbulence  model  to  account  for  the  enhancment  of  vertical  mixing  by  LC.  This  included 
modifications  to  the  representation  of  the  vertical  turbulent  momentum  fluxes,  which  are  used  in 
the  mean  horizontal  momentum  equations  as  well  as  in  the  turbulence  model’s  TKE  and  TLS 
equations,  and  to  the  stability  functions  of  the  turbulence  model. 

Both  the  KC04  and  H15  LCMPs  were  implemented  in  the  MYL2.5  turbulence  model  used  in 
NCOM  and  tested  by  conducting  simulations  of  the  simple,  light-wind,  mixing  case  used  by  MW97 
(and  also  by  KC04  for  their  own  testing),  the  upper-ocean  thermal  structure  at  OWS  Papa  during 
1961,  and  Hurricane  Ivan  in  the  Gulf  of  Mexico  in  2004. 

The  simple  test  case  used  by  MW97  and  KC04  consisted  of  forcing  by  a  constant  surface  wind 
stress  of  0.037  Pa  (which  corresponds  to  a  wind  speed  of  about  5  m/s),  a  small,  cooling  surface  heat 
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Fig.  20  —  Ratio  of  maximum  values  of  Km  in  SML  for  Hurricane  Ivan  simulations  with  the  H15  LCMP  and  without 
a  LCMP  at  09Z  15  September  2004.  The  white  vectors  show  the  surface  winds. 
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Fig.  21  —  Ratio  of  maximum  values  of  Km  in  SML  for  Hurricane  Ivan  simulations  with  the  KC04  LCMP  and  without 
a  LCMP  at  09Z  15  September  2004.  The  white  vectors  show  the  surface  winds. 
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flux  of  5  W/m2,  and  an  initial  MLD  of  33  m.  The  use  of  the  KC04  LCMP  increased  the  maximum 
rate  of  mixing  in  the  SML  by  about  a  factor  of  two  from  about  200  to  400  m2/s,  and  increased  the 
MLD  by  about  3  m.  The  increase  in  the  maximum  rate  of  mixing  from  LC  was  roughly  consistent 
with  the  results  in  MW97  and  KC04.  The  use  of  the  H15  LCMP  resulted  in  a  similar  increase  in 
the  MLD  of  about  3  m,  but  increased  the  maximum  rate  of  mixing  in  the  SML  to  a  value  of  about 
1500  m2/s. 

For  the  Papa  simulations  for  1961,  no  direct  wave  observations  were  available.  Hence,  the  wave 
held  was  estimated  from  the  Papa  winds  by  establishing  a  relation  between  the  local  wind  speed 
and  the  wave  spectrum  based  on  9  years  of  wind  and  wave  observations  from  NO  A  A  Buoy  46005, 
located  in  the  northeast  Pacific  about  1100  km  SE  of  OWS  Papa. 

The  use  of  the  KC04  LCMP  at  Papa  increased  the  rate  of  mixing  within  the  SML  in  the  spring 
and  summer  by  about  a  factor  of  2  over  the  original  MYL2.5  MLM  without  a  LCMP.  The  predicted 
MLD  at  Papa  was  slightly  increased,  but  not  enough  to  provide  good  agreement  with  the  observed 
MLD  in  the  summer  and  fall  of  1961.  However,  addition  of  the  L94  parameterization  of  unresolved 
mixing  processes  (Large  et  al.  1994)  deepened  the  ML  enough  to  provide  good  agreement  with  the 
observed  MLD  at  Papa.  This  was  consistent  with  KC94,  who  reported  that  their  implementation 
of  L94  in  the  MYL2.5  turbulence  model  significantly  improved  the  predicted  MLD  at  Papa. 

The  simulation  at  Papa  with  the  H15  LCMP  resulted  in  a  fairly  good  prediction  of  the  observed 
SST,  MLD,  and  ISODs  without  the  need  to  incorporate  additional  mixing  mechanisms,  e.g.,  such 
as  the  L94  parameterization.  The  H15  LCMP  increased  the  rate  of  mixing  within  the  SML  in  the 
spring  and  summer  at  Papa  by  about  a  factor  of  2-5  over  the  rates  of  mixing  predicted  by  the 
KC04  LCMP;  hence,  about  a  factor  of  4-10  over  the  original  MYL2.5  MLM. 

For  the  simulations  of  Hurricane  Ivan  in  the  Gulf  of  Mexico,  the  use  of  the  KC04  and  H15 
LCMPs  generally  resulted  in  significantly  increased  mixing  rates  in  the  SML  near  the  hurricane, 
i.e.,  the  mixing  rates  are  increased  by  a  factor  of  2  to  3  for  KC04  and  by  a  factor  of  5  to  10  for  H15. 
The  increased  mixing  rates  result  in  a  significant  reduction  in  the  vertical  shear  of  the  wind-driven 
currents  in  the  SML,  which  results  in  a  reduction  of  the  surface  value  of  the  wind-driven  current 
near  the  hurricane  of  20  cm/s  or  more.  The  depth  of  mixing  is  also  generally  increased,  which 
results  in  increased  cooling  of  the  SST.  For  the  KC04  LCMP  the  additional  SST  cooling  is  small, 
i.e.,  generally  less  than  0.2°C;  but  for  H15  the  increased  cooling  is  greater,  i.e.,  0.2-0.4°C. 
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CALCULATION  OF  SDC  FOR  A  MONOCROMATIC  WAVE 


The  Stokes  drift  current  (SDC)  Vs  can  be  computed  for  a  monocromatic  (i.e.,  single  frequency) 
wave  using 

Vs(z)  = 

where  a  is  the  wave  amplitude,  k 
(tanh(kH)g/k)2  is  the  phase  speed, 
and  z  is  the  distance  to  the  surface. 

In  deep  water,  this  simplifies  to 

Va( 

where  the  phase  speed  is  cp  =  (g/k)  2 

For  a  wave  field  consisting  of  waves  of  different  amplitudes,  frequencies,  and  directions,  the 
contributions  to  the  SDC  from  the  different  waves  can  be  summed  (in  a  vector  sense). 

The  wave  amplitude  a  is  related  to  the  wave  energy  E  as 

E  =  ia2.  (A3) 

Note  that  the  wave  energy  given  by  Eq.  (A3)  has  units  of  length  squared,  which  is  a  common 
practice  when  discussing  surface  waves.  Multiplication  by  the  water  density  and  the  acceleration 
of  gravity  gives  the  proper  units  of  energy  per  m2  of  surface  area. 


2  cosh(2  k(H  +  z)) 

-  ( akYCn - TT, - : — 

2  skill2  (/cR) 

=  27r/L  is  the  wavenumber,  L  is  the  wavelength,  cp  = 
g  is  the  acceleration  of  gravity,  H  is  the  total  water  depth, 


z)  =  (ak)2cp  exp(2/cz),  (A2) 
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CALCULATION  OF  SDC  USING  WAVE  SPECTRA  FROM  NO  A  A  BUOY 


Some  of  the  NOAA  buoys  provide  observations  of  both  the  surface  wind  velocity  and  the  wave 
energy  spectra.  These  data  allow  the  calculation  of  the  average  wave  spectra  at  different  wind 
speeds.  The  average  wave  spectra  at  the  different  wind  speeds  can  then  be  used  to  estimate  the 
SDC  for  a  given  wind  speed. 

The  spectral  wind  and  wave  data  that  were  used  here  were  from  NOAA  Buoy  46005,  which  is 
located  in  the  northeast  Pacific  at  131.001°W,  46.100°N,  which  is  about  1100  km  east-southeast 
of  the  location  of  OWS  Papa.  Since  the  location  of  the  buoy  is  somewhat  similar  to  that  of  OWS 
Papa,  the  wind  and  wave  conditions  should  be  roughly  similar. 

Data  from  Buoy  46005  for  the  years  1996  through  2004  were  used.  For  these  years,  the  spectral 
wave  data  from  the  buoy  consist  of  38  values  of  the  wave  energy  at  frequencies  from  0.03  to  0.4  Hz 
at  0.01  Hz  intervals.  Note  that  no  directional  information  is  provided  for  the  wave  data  from  Buoy 
46005.  Averages  of  the  spectral  data  were  computed  within  1-m/s-wide  wind-speed  bins  from  1 
to  20  m/s.  No  winds  above  20  m/s  were  recorded  at  Buoy  46005  during  the  9- year  period.  The 
wind  speed  used  to  characterize  the  wave  spectrum  in  this  calculation  is  the  wind  speed  at  the  time 
of  the  wave  observation.  An  alternative  would  be  to  use  some  time-weighted  wind  speed  over  the 
previous  few  hours,  but  this  was  not  done. 

Figure  B1  shows  the  averaged  wave  spectra  from  NOAA  Buoy  46005  for  wind  speeds  of  1,  5, 
10,  15,  and  20  m/s.  The  wave  energy  increases  with  increasing  wind  speed.  There  are,  however,  a 
couple  of  notable  aspects  to  the  wave  spectra  in  Fig.  Bl:  (1)  the  location  of  the  peak  of  the  wave 
energy  spectra  shows  little  change  with  wind  speed,  and  (2)  there  is  a  notable  amount  of  wave 
energy  present  at  low  wind  speeds  of  1  to  5  m/s. 

There  are  several  limitations  to  estimating  the  SDC  from  a  wave  spectra  based  on  the  current 
local  wind  velocity.  The  wind  duration  and  fetch  and  the  propagation  of  waves  from  other  areas 
cannot  be  accounted  for.  Additionally,  variation  in  the  direction  of  the  waves  making  up  the  energy 
spectrum  cannot  be  accounted  for;  hence,  the  direction  of  the  waves  and  the  SDC  is  assumed  to 
be  downwind. 

The  SDC  profiles  estimated  from  the  wave  spectra  at  Buoy  46005  for  different  wind  speeds  were 
compared  with  the  SDC  profiles  computed  from  the  observed  wave  spectra  using  data  for  the  year 
2000,  and  the  mean  and  rms  differences  versus  depth  are  shown  in  Fig.  B2.  The  SDC  estimated 
by  this  method  will  be  referred  to  in  this  report  as  the  Buoy  Wave  Spectra  (BWS)  SDC.  Plots  for 
the  same  errors  computed  using  data  from  buoy  46005  for  other  years  between  1996  and  2004  (not 
shown)  look  similar  to  Fig.B2. 
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Wave  Energy 


Fig.  B1  —  Averaged  wave  spectra  from  NOAA  Buoy  46005  for  wind  speeds  of  1,  5,  10,  15,  and  20  m/s.  The  wave 

energy  increases  with  wind  speed. 
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Appendix  C 


PARAMETERIZATION  OF  UNRESOLVED  MIXING  PROCESSES  BY  LARGE 

ET  AL.  (1994) 


As  noted  in  the  introduction  to  this  report,  tests  of  upper-ocean  mixing  in  the  open  ocean 
with  turbulence  models  such  as  the  MYL2.5  MLM  have  frequently  found  that  the  models  do  not 
predict  as  much  mixing  as  is  observed.  It  has  been  hypothesized  that  is  due,  not  so  much  to  the 
fact  that  the  mixing  mechanisms  in  the  models  are  incorrect,  but  that  some  mixing  processes  are 
either  unrepresented  or  under-represented,  including  sources  of  background  shear  (e.g.,  internal 
waves  and  inertial  gravity  wave  pumping),  surface  waves,  and  Langmuir  circulation. 

To  account  for  shear  instability  mixing  below  the  SML,  Large  et  al.  (1994,  referred  to  in  this 
report  as  L94)  proposed  extending  such  mixing  to  Richardson  numbers  Ri  above  the  normal  critical 
value  of  0.2  to  0.3,  where  shear  instability  is  suppressed  by  the  strength  of  the  density  stratification. 
The  L94  mixing  enhancement  extends  the  mixing  to  Ri  =  0.7  and  is  described  by 

K0  Ri  <0 

Kl9 4=  Ka(l  —  (Rj/0.7)2)3  0  <  Ri  <0.7 

0  Ri  >  0.7, 

where  Ka  =  50  cm2/s.  In  NCOM,  XL94  is  computed  from  Ri  and  is  applied  to  the  vertical  mixing 
of  the  momentum  and  scalar  fields  (e.g.,  temperature  and  salinity)  by  using  the  maximum  of  the 
usual  vertical  mixing  coeffient  {Km  or  Ku  for  momentum  or  scalar  fields,  respectively)  and  Kl 94 
to  compute  the  vertical  mixing. 

This  scheme  was  utilized  by  L94  in  conjunction  with  an  adaptation  of  the  atmospheric  bound¬ 
ary  layer  model  of  Troen  and  Mahrt  (1986)  to  the  ocean  (known  as  the  KPP  MLM),  and  by  Kantha 
and  Clayson  (1994,  referred  to  in  this  report  as  KC94)  in  conjunction  with  the  MYL2.5  turbulence 
closure  model.  Both  L94  and  KC94  found  that  the  addition  of  this  mixing  improved  the  agreement 
of  predictions  of  the  ocean  surface  mixed  layer  with  observations  from  several  open-ocean  data 
sets,  including  data  from  OWS  Papa.  However,  a  question  is  whether  such  a  mixing  enhance¬ 
ment  is  needed  in  shallow,  coastal  water,  where  overmixing  is  sometimes  more  of  a  problem  than 
undermixing. 
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